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ABSTRACT 


The analytical procedure for computing second-order perturbations in 
rectangular coordinates, according to Brouwer's theory of planetary motion, 
is given. Single- and double -harmonic analyses and the multiplication of 
Fourier series with numerical coefficients are used in the computations. In 
the series multiplication, a variable tolerance is considered, enabling us to 
avoid the difficulties arising from a small divisor. 

Also presented is an example computing that part of the second-order 
perturbation of Mars containing the masses of Jupiter and Saturn. The 
analytical solution of this perturbation is compared with the numerical inte- 
gration of the differential equations defining this perturbation. The 

numerical integration covered the interval from 0 to 40, 000 days. The com- 

-9 

parison shows an agreement within 1 X 10 

RESUME 

La procedure analytique de calcul des perturbations de second ordre en 
coordonnees rectangulaires, d'apres la theorie de Brouwer des mouvements 
planetaires, est exposee. Des analyses d'harmoniques simples et doubles et 
la multiplication de series de Fourier avec des coefficients numeriques sont 
utilisees dans les calculs. Dans la multiplication de series on tient compte 
d‘une tolerance variable qui permet d'eviter les difficultes provenant d'un 
petit diviseur. 

Un exemple de calcul de la perturbation de second ordre de Mars due aux 

masses de Jupiter et Saturne est egalement presente. La solution analytique 

de cette perturbation est comparee avec 1 ! integration numerique des equations 

differentielles qui definis sent cette perturbation. L’ integration numerique 

couvre l f intervalle de 0 a 40. 000 jours. La comparaison montre une 

-.9 

concordance de 1 X 10 


v 



Pesjone 


MajiaraeTCH aHansiTH^iecKHM cnocoO BbiHMCJieHHH BoaOyacfleHHH 
BToporffl nopiiflKa b npnMoyrojiBHtix KOopflHBBTax , coxviacHO Bpayap- 
ckoh TeopHH miaHeTHoro ABHxeHsiHo 0#ho- h flByx-rapMOHnaecKHe 
aHajiH3bi h yMHOxeHne ceprni $oypi>epa c hh c ji o bmmh K03$$HpHeHTaMH 
HcnojitaoBaHH b pacaeTax . ITpH yMHoace hum cepiiH Obuia yHHTaHa ne- 
peMeHHan TOJiepaHpuH, iio3bojihh nafiexaTB Tpy#HO cth , B03HHKaro- 
mvie SjiaroflapH MajiOMy paa^eaHTeaio. 

TaicEe npuBeaen npHMep BbrencjieHHa MapcoBux B036yxja,ewivi bto- 
poro nopn^Ka, coflepxamiix Maccw ffinwTepa h CaTypHa • AHajiHTHaec- 
Koe peiueirae othx B 036 ys:fleHHH cpaBHeHO c ancjiOBoft HHTerpauHew 
AH$$epeHiiHajiBHbix ypaBHeHHH, onpeflejiaiomHx .3 tux B036yx,zi;eHHft . Hhc- 
jiOBaa HHTerpautaa noicpbuia HHTepBfiui ot 0 ao 40,000 AHeii. CpaBHeHHe 
noKa 3 biBaeT corjiacae b npeaejxax I x 10 


VI 



SECOND-ORDER PLANETARY THEORY 


Part I: Outline of the Method 

S. E. Hamid 

L INTRODUCTION 

The author has successfully applied Brouwer's theory of general pertur- 
bation in rectangular coordinates to obtain a first-order planetary theory for 
all the principal planets except Pluto (Hamid, 1968). The advantage of 
Brouwer's theory over other planetary theories is its convenience when higher 
order perturbations are considered. 

In this report, the adaptation of the theory in the computation of second- 
order perturbations is discussed. General computer programs have been 
developed for the computation of the different second-order terms of plane- 
tary perturbations. These programs have been applied for the planet Mars 
to compute the second-order perturbations factored by the product of the 
masses of Jupiter and Saturn. The numerical results obtained have been 
tested successfully against the numerical integration of the differential equa- 
tions satisfying these perturbations. 


This work was supported in part by grant NGR 09-015-002 from the 
National Aeronautics and Space Administration. 
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2. THE EQUATIONS OF MOTION 


Consider a set of rectangular axes, the x axis corresponding to the 
direction from the sun to the perihelion of the orbit of the perturbed planet 
at a given epoch, and the z axis perpendicular to that orbital plane at this 
epoch. Then, the perturbations 6x, 6y, 6z in the rectangular coordinates 
satisfy the following set of differential equations: 


d 2 Sx Sx X 0 , c c » „ 

+ P —r 7 — (x Q Sx + y Q 6y) = G 

dt r 0 r 0 


X 


d 2 Sy 6y 3|ay 0 , 

— 2 + I 1 3 " 5 ( x 0 6x + Yq 6 y) = G 

dt r o r o y 


d 6z , 6z ^ 

~T~ + i" T = G z • 


dt 


0 


( 1 ) 


The quantities (x^, y^, z^) are the coordinates of the planet, with its unper- 
turbed orbit assumed at epoch, while r^ denotes the heliocentric distance 
of the planet. 


The functions G^, G^, can be separated into different parts of 
descending order of magnitude, the first part giving rise to first-order per- 
turbations, the second part to second-order perturbations, and so on. We 

and let 6x, , 6'y, , 


denote these parts by G, , G, , G, ; G 0 
r 7 lx ly lz 2x 


6z^ ; 6 x 2 , 5y 2 , 6z 2 ; 
coordinates . 


JX^ 0 ,G 0 , G ) • • • « tfeiJ-V-L j ^ y 1 

ly lz 2.x 2y 2z 1 7 1 

be the corresponding perturbations in the rectangular 


The first-order perturbations 6x^ , 6y^ , 6z^ will satisfy equations (1) 

when the values of G , G , G are put equal to G, , G, , G , , and 

x y z lx l y 1 z 

similarly for higher order perturbations. 
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We consider the second-order perturbations, which are written as 


d 6 x 2 3jiXq 

”2 + — - — 5 ~ (x 0 6x 2 + ^0 6 Y2 > = G 2x ’ 


dt 


r o r o 


d 6y ? 6y 3|iy 0 

2 + M- — ~ — 5— ( x 0 6x 2 + YO d Tz) = G 2y , 


dt 


r r 

0 0 


d 6 z_ 6z_ 

— — 2 + fJ. 3“ = G 2 z • 

dt r Q 


( 2 ) 


The solution of equations (2), given by Brouwer and Clemence (1961), takes 
the following form: 
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(3) 


6z_ = q. 


J 


G 0 dt 
2z 


q 


/ 


^2 G 2z dt 


For the definition of the different partial derivatives 
Yq and of the quantities q^, q^, in equations (3), see 
( 1961 ). 


of the coordinates x^, 
Brouwer and Clemence 
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3. THE DEVELOPMENT OF G 0 , G. , G 0 

Zx Zy Zz 

The expressions for G ^ > G^ ^ ave following form: 


a 2 R r 


a 2 R 


G 


o_ 

k^k 


a 2 R, 


2x ^2 9 X k ^ 3x u 3y u 9y k ^ Bx^. 3z t ~“k 


5z, 


k k 


n /8 2 R n 9 2 R n 9 2 R \ 

2 ^8x k a Xj 6x j + 8x k 8 Yj 6y j + 8x k 8z . 6z j j 


+ M- 


9 x k 1 5\ K 2 / Jk . X k y k 

2 T -TT i 6x k + 3— - 15 — 
r k r k/ V r k r t 


6 x k 6 yk 


+ | 1 fk.BVt k 2 + 3i 6z 2 

+ 1 2 5 2 7 I y k + 2 5 0 k 

r k r k / r k 


2 

a r. 


6x. + 


* 2 *o 6 /»» 6 

Oy. + - — - — oz, 


2y 9y k ax k -k 8 2 'k 3y k 8z k ~k 


n / rv2 


. » 2 R ° ■ 
— y 8 y k ax j X j ay k 9y j y j + 8 y k 9z j 


+ M- 


ilk 15 Vk 
2 5 " 2 7 

r k r k 


6x, + 3-4-15 


X k y k 


6x k ^ y k 
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. 9 15 \ r 2 3 r 2 

+ 1 2 T " T — 6y k + 2 “5 6z k 
r k r k / r k 


. 82r o , 1)2 R o . , 82r o 

t 2z 3z 3 x x k + 3z, 3y ^k + A 2 k 

k k k 7 k 3z, 

k 


g2 R o 1)2 Rp 

^L> \ Sz k 8x j 6x j + 8z k 8y j 6y j + 3z k 3z i 6z j 




j=l 


1 Z . J * 

k J 


+ t* l 3 4f 6x k 6z k + 3 % 6y k 6z k I • 

r k r k 


(4) 


In equations (4), we have 


x. 

J 


z . 
J 


= the rectangular coordinates of the disturbing planet j, 
with its unperturbed orbit assumed at epoch; 


x. 


k’ y k’ Z k 


6x k , 6y k , 


= the rectangular coordinates of the disturbed planet k, 
with its unperturbed orbit assumed at epoch. Note that 

x k 1 V y k * y o- z k = 0i 

6z^ = the perturbations in the rectangular coordinates of 
the disturbed planet k; 


6 Xj, 6y^, “ the Perturbations in the rectangular coordinates of 

the disturbing planet j; 

2 

p = k (1 H-m^b where k is the gaussian constant and is the mass 
of the disturbed planet k; 

Rq = the well-known disturbing function of the different disturbing 
planets on planet k, given by 
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(5) 



where is the mutual distance of planets k and j, and r. is 
the heliocentric distance of planet j. 


We note that 23 represents the sum over all the disturbing planets j. For 

j^k 

example, if we consider the theory of Mars, then we have k = 4, and j will 
take the numbers 1, 2, 3, 5, 6, 7, 8, corresponding to the effects of Mercury, 
Venus, Earth, Jupiter, Saturn, Uranus, and Neptune. In this report, we shall 
exclude the effect of Pluto. 


For the perturbations 6x, , 6y. , 8z, and 6x., 6y., 6z., we shall consider 

k k k j j j 

the values derived from the first-order theory. The Sx^., Sy^, 6z^ are com- 
posed of different parts, owing to the perturbations of the different disturbing 
planets . 


If we let Fj k (jk), F^^jk), F 3 ^(j k ) be, respectively, the first-order per- 
turbations in 6x^, 6y k , 6z^ due to the disturbing planet j, then, 

6x k = 2 F ik (jk) • 

k* 

6 ^k = Z F 2k (jk ) » 

kk 

6z k = Z F sk ( k) • (6) 

kk 

Similarly, 

6x j - S v ij) < 

i* j 
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(T) 


6y j = E F 2j (i « • 

6z. - £ *««> ' 


The coefficients of Sx^., in equations (4) can be written as follows: 


3 2 R 


2 

a &, 


ax, 


12 ' t ’ x i° k) • 8y^ = 12 °yi <ik) 




2 

o R, 


8x k 8y k 


j*k 


* 2 4 ’s2 (ik) ’ 


a 2 R 0 ^ 


a 2 r. 


8z k 8x k 


82r q 

9z k 9y k 


= 2 **i (Jk) * 


j*k 


= ^^2^^ * 


j^k 


82r 0 

8x k 8z k 




j*k 


where 


2 

d R r 


8y k 8z k 1 


= 12 *y3 


j^k 


a 2 r, 


3z 


= S *z3 (jk) ’ 


j* k 


(8) 


+xl (jk) = k m. 


+ 


k j 


3(x i - x k > 

<■ 

k J 


4> x2 (J k ) = k m j 


2 3 |x i ~ x k> (y i ~ y k’ 


A k‘ 

k J 


*x3 (jk) = 


k m. 3 (x. - x, )z . 

, i j k j . 


a 3 
k j 


c), y2 (jk) =k m. 


i . 3(y i - y k } 


l A kj 


aI 

kj 
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4 > 

4> 

4> 

<t> 


z3 (jk) =k 2 m. 

yi (j k ) = * x2 (jk) 

zl (jk) = *x3 (jk) 
z2 (J k ) = 4> 3 (jk) 



i 


(9) 


The coefficients of 6x., 6y., 5z. in equations (4) can be written 

J J J 


8 2 R 0 

3x,3x. 9 xl ( j k) * 
k J 


3 2 R„ 


9 y k 9x i ~ e y 1 ^ 


a 2 R, 


art: = 9 ziW 

k J 


2 

3 R, 


3x, dy. x2 
k 


= ®„o(jk) , 


2 

3*R 


Kr=e,W . 


9y k 9y j 72 


8 2 R 0 

9z. 9y. " 6 z2^ k * ’ 

k 'j 


3 2 R, 


3x^3zj x3 


where 


= e „^(jk) , 


2 

p 

9y^j = 0 y3 (jk) 


e x l(jk) = k m. 


kj 


_ 1 _ 

3 
r . 

J. 


, 2 

3x. 


5 
r . 
J 


8 2r 0 „ 

9z,9z. 0 z3^ k * • 

k J 


(10) 


3(x - Xj/ 

4“ 

kj 


0 x2 ( jk) = k m. 


3(x. - x k )(y. - y k ) 3x. y . 


kj 
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9 x3<S k > = k m i 


3(x. - x, )z . 3x. z . 

1 k I + i. ,J 


kj 


5 
r , 
J 


e y2 (j k ) = k m. 


1 _1_ 

.3 " 3 

A.. r. 
k J J ■ 


, 2 ' 
3 y^ 


3(y j - ^ - + _ 

A?. r 5 

k J J 


0 y3 (jk) = k m.. 


3 (y, - yjz: 3y z. 

- - -- j + — L 


kj 


5 
r . 

J 


0 z3 (jk) = k m. 


1 J_ 

l A, . r . 

k J J 


3z[ 

3z 2 

* 5 

*• s~ 

Af , 

r . 

kj 

j - 


9 yl (jk) = 0 x2 (j 


9 zl (jk) = e x3 (jk) * 


e z2 (jk) = 0 y3 (jk) 


(ID 


2 2 2 

Finally, the coefficients of (6x^_) , (6y^) , (6 z^) , 6x^ Sy^, Sx^-Sz , 

5y 6 z v in equations (4) can be rewritten as follows: 

JbC 


4* xl (k) - H 4.5 -1 - 


3 

7.5 -j 

r k 


x k y k 


^x2( k ) =M 3 -f - 15 -T 
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*k 


x k y k 


^(kJ^U.S-f- 7. 5^ 


^ 4 (k)= M- 5 ^ 


^ yl( k ) " ^ l 1 ’ 5 5 

X r k 


7. 5 


2 

Vk' 


x, 


2 

Vk 


V k,a# “l 3 "f ■ 15 ^ 


4, y3 (k)= k(4.5^| - 


3 

y k 
7 , 5 -j 

r k 


4J y4 (k) = M- 5 — 


x r 


+ zl (k)= ^ 3 ~ 


^ Z 2 (k) = ^ 3 ^t 


(12) 
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With the above definitions of the various coefficients in equations (4), we 
have the following: 


G 2 X = E »*,<*> + E 4 X?- C kl E F 2k <J ' k) 

j^k j#k j^k j#k 


+ E + X 3« k > E F 3k« k| 

j*k j*k 


+ E 


xl 


« k, E 


F :;( ii » 1 ° x2 


E F 2j (ij) + 6 x3« k) E F 


3j 


,(ij) 


i*j 


i*j 


+ + xl (k) 

E F lk« k > 

_j*k 

2 

+ 4 x2 (k) 

E F ik« k > 

_>k 


E F 2k<i k > 

J* k 

+ + X 3< k > 

E F 2k < i k) 

J* k 

2 

+ + x4 (k) 

E F 3k<i k) 

j*k 

2 

4 


G 2y ' E ♦yl (ik) E F lk (jk) + E *y2 (Jk) E F 2k< 

j^k j*k j#k j^k 

+ Ev jk) E F 3k« k) 

j^k j^k 


+ E 

j#k 


yi 


(iwE 


F lj«> +e y 2 


<j k) E 


F *J«> 


+ 0 


y3 ,jk) E 


F 3j"J> 
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~ “ 

2 

- -1 


r- - 

+ K yl (k) 

2>ik (jk) 

+ + y2 (k) 

L r ik<i k > 


E F 2k'i k » 


J*k 


Lj*k 


J* k 

+ V k> 

E F zk« k) 

2 

+ V k > 

E F 3k<i k > 

2 

3 


Lj* k 


J*k 



g 2z = E Ki <*> E F ik« k > + E E F 2 k « k) 

j^k j^k j^k j*k 


j*k j^k 


♦ei 

j^k 


a ( *>E 

i*j 


F 2 .(ij) + 0 


z2 


<«E 

i*j 


F 2jW> + 0 z3 


e 

i*j 


F 3j (ij) 


+ +z!< k > 

E F ik« k > 


E F 3k<J k > 


j^k 


j* k 


+ *z2< k > 

E F 2k« k > 


E F 3k« k > 


J* k 


_j*k 


Let us now look more closely at the different terms defining the quan- 
tities G 9 , G 9 , G 9 . The c|> and 0 terms can be represented by double 
eh x cy 

Fourier series in the mean anomalies f i ^ of the disturbed planet k and 

the disturbing planet j. These series can be obtained by computing the 
special numerical values of cf> and 0 for different combinations of equidistant 
values of the mean anomalies i^ 9 4 y These special. values are then subjected 
to double -harmonic analysis. 
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The 4* terms can be represented by Fourier series in one argument, 
the mean anomaly f ^ the disturbed planet k. These series can be 
obtained by computing the special numerical values of ij j for different 
equidistant values of the mean anomaly i ^ and then subjecting these 
values to single -harmonic analysis. 

In other words, by expressing the c|>, 0, and 4 J terms as Fourier series 
in the mean anomalies, we can avoid analytical expansions. Only double- 
and single -harmonic- analysis techniques can be applied. This is what we 
have done in the present work. In fact, a general computer program can be 
constructed to have as output the Fourier representations of the different 
cj>, 9, and 4* terms for any given values of j, k. 


The terms F^(jk), F^k^l^ have already been obtained in the 

first-order theory. It should be remembered that these perturbations in 
rectangular coordinates are composed of two parts: the periodic and the 
secular. The periodic part is represented as double Fourier series in the 
mean anomalies i y 1 and the secular part by the product of the time t 
(measured from the given epoch) and a single Fourier series in the mean 
anomaly Let the periodic part be denoted by f^(jk), ^k^^* f^^(jk), 

and the secular part, by t S^^(k), t t Hence, 


F ik^ k ^ = f ik^ k * + t 5 



(14) 


where i = 1 ,2, 3. 


Let us consider the part 53 <|> 


K xl 


From the above remarks, this part 


is represented by the summation of different Fourier series, and each series 

is represented in the mean anomalies l . and i, . For example, if we are 

J ■ y' 

considering the theory of Mars, we have k = 4, and Z-r <(> , (jk) will be eom- 

j^k xl 

posed of the sum of seven Fourier series: the first series ini^, the 

mean anomalies of Mercury and Mars; the second series ini^* ^ the 

mean anomalies of Venus and Mars; and so on. Similarly, the part^ F , (jk) 

j*k IK 

is composed of the sum of different Fourier series, and each series is 
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expanded in the mean anomalies i y In addition to these Fourier series, 

this part contains a term t multiplied by a Fourier series in single argument 
i k* the mean anomaly of the disturbed planet* In fact, for i = 1, 2, 3, 


i>ik<*>=L f ik« k)+ts ik< k > . 

j#k j^k 

where 


(15) 


s ik< k > - E v k) • 

j*k 

We note that t S^(k), t S^(k) 9 t S^^(k) are the total secular perturbations in 
rectangular coordinates of all the disturbing planets on planet k. 

Similar considerations apply for the different parts of equations (13). 
Hence, G^, can be represented by the following equations: 


G 2x= G 2x ktt< k > ' ‘ 2 + E [G 2 X jkt« k > ' t+ G 2*jk« k >l 

j*k 


+ G 2xmjk (m j k) ■ 

j m 


°2y G 2yktt^ ' t + 5^ G 2yjkt^ ‘ 1 + G 2yjk^ k ^ 

j*k 


+ EE C 2ymjk ,mjk) ' 

j m 
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( 16 ) 


G 2z = G 2zktt< k > ' ‘ 2 + Z [G 2zjkt<i k > ' ‘ + G 2zjk ( j k) l 

j*k 


EE 


G 2zmjk (m£k) 


m 


where G(k), G(jk), and G(mjk), appearing on the right-hand side of equations 

(16), denote, respectively, Fourier series in one argument, the mean 

anomaly i , ; in two arguments, the mean anomalies i ., i, ; and in three 

^ j v' 

arguments, the mean anomalies i , i ., i. . The double summation L L 

m J ^ j m 

means that m and j take the values corresponding to all the disturbing 

planets, excluding m = j and avoiding double counting. 


Following are the expressions for the different G(mjk) in equations (16): 

3 

G 2qmjk (mjk) = E [<,> qi (j;k) f ik (mk) + V" 1 * 1 £ ik (jk) + 0 qi (jk) £ ij (mj) 

i= 1 

+ e qi (mk) f im (jm)] + 2^ ql (k) f lk (jk) f lk (mk) 

+ ^ q2 (k)[f lk (jk) f 2k (mk) + f Ik (mk) f 2k (jk)] 

+ 2+ q3 (k) f 2k (jk) f 2k (mk) + 2+ q4 (k) f 3k (jk) f 3k (mk) , 

where q = x, y, and 

3 

G 2zmjk (m£k) “E^zi 1 ^ f ik (mk) + + zi (mk) £ ik<j k) 
i= 1 

+ 0 . ( jk) f..(mj) + 0 .(mk) f. (jm)] 
zi J ij J zr im 
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+ + z i(k)[ f i k (j k ) f 3k (mk) + f lk (mk) f 3k (jk)] 

+ ^ z2 (k)[f 2k (jk) f 3k (mk) + f 2k (mk) f 3k (jk)] . (17) 

The terms G(jk) (not multiplied by t) in equations (16) are expressed as 
follows: 

3 

G 2qjk (jk) = S [4> qi {jk) £ ik {jk) + 9 qi ( j k) f ij (k ^ + * ql (k) f lk (jk) 
i= 1 

+ ^q2^ f lk^ k ^ f 2k^ k ^ + ^ q3^ f 2k^ + ^q4^ f 3k^ ’ 

where q = y, and 

3 

°2zjk (jk) £ ik (jk) + 9 zi (ik) f ij (kj)1 + + z 1 k £ lk (jk) f 3k<j k) 

1=1 

+ 4> z2 ( k ) f 2 k (j k ) f 3k (j k ) • ( 1 8) 

The terms G(jk) (multiplied by t) take the following forms: 

3 

G 2qjkt<« = E B-qifW S ik (k > + S ij<3>] + 2 + q l< k > £ lk« k > S lk< k > 

i=l 

+ '*' q2 ' k « £ lk« k) S 2k< k > + £ 2k ( j k) S lk< k » 

+ 2 ^ q 3^ f 2k^ S 2k^ + 2 ^q4^ f 3k^ S 3k^ ' 
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where q = x, y, and 


3 

G 2zjkt<i k » = 2Xi<*> s ik(W + e zi(j k > 

i= 1 

+ + 2l (k)[f lk (jk) s 3k (k) + f 3k (jk) s lk (k)] 

+ ^ z2 (k)[£ 2k (jk) S 3k (k) + f 3k (jk) S 2k (k(] (19) 

For the terms with coefficient t in equations (16), we have 

G 2 q kt,< k > = + q l< k > S ?k (k > + %2 (k > S lk< k > S 2k< k > 

+ + q 3< k)S 2k <k) + + q 4 (k,S !k' k > • 

where q = x, y, and 

G 2zktt< k > = + zl < k > S lk< k) S 3k< k > + S 2k< k) S 3k< k > ■ < 20 > 


In equations (17) to (20'), all the different terms on the right-hand side 

are expressed in Fourier series in the mean anomalies i , i i, . We 
■ m y k 

have already outlined how we obtain these series. We are now in a position 
to evaluate the Fourier representations in mean anomalies of the functions 

G 2qmjk {mjk) ’ G 2qjk (jk) ' G 2qjkt’ G 2qktt’ where q denotes the parameters 
x, y, and z. 


Let us consider, for example, G^q^^mjk), whose expressions are 
given in equations (17). Since we have the Fourier series for all the terms 
appearing in the right-hand side of equations (17), we can, by the technique of 
multiplying Fourier series, obtain the Fourier series representing 
G 2 qmjk^ m j k ^ case » we did not resort to triple-harmonic analysis 

because it would have been excessively laborious. In fact, we constructed 
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a general computer program that has as input the numerical values of m, j, 
and k and that will give as output the Fourier representations of ^qmjk^'M 
in the mean anomalies & y 4^ (q denotes the values x, y, and z). 

In computing the Fourier representations of C^qjktO^ and C^qjkUk) 
for q = x, y, z, we can use the double -harmonic- analysis approach or the 
multiplication-of-series approach. To compute the Fourier representations 
of G^qktr^) 4° r P = y> z » we can very conveniently use the single- 
harmonic-analysis technique. 


In our work, we have a general program that, for given j, k as input, 
produces as intermediate output the Fourier representations of 


G 2qjk (jk) ’ and G 2qktt (k) f ° r * = Z - 
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4. THE DECOMPOSITION OF 6x 2 , Sy 2 , 5z 2 

In the previous section, we developed the different Components of the 
functions G 2x , G^, G,,^. We found that these functions are generally com- 
posed of the summation of the following series; 


A. Fourier series in three arguments. 

B. Fourier series in two arguments. 

C. Fourier series in two arguments, multiplied by the time t. 

2 

D. Fourier series in one argument, multiplied by t . 


By substituting the general expressions of G 2x> G^, G 2z in equations 
(3), we can see that 6 x2 , ^y2* ^z2 ^e composed of the following 

different parts, where q takes the values x, y, z: 


A. 

B. 

C. 

D. 

E. 

F. 

G. 

H. 


Fourier series in three argument (i , f •> denoted by 

Xxx j 


Fourier series in two arguments (i y l ^), denoted by 6q 2 j^.(jk). 

Fourier series in two arguments (f ., ) multiplied by the time 

J K 

t, denoted by 6 q 2 j kt (jk). 

Fourier series in one argument (i k ) denoted by 5q2^(k). 

Fourier series in one argument (1^) multiplied by t, denoted 

by 6q 2kt (k). 

Fourier series in one argument (i^) multiplied by t , denoted by 

5q 2ktt 3 

Fourier series in one argument (i^.) multiplied by t , denoted 

b ^ 5< !2kttt (k) - 4 

Fourier series in one argument (l k ) multiplied by t , denoted 

b y ^2kt t tt (k) - 
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We must remember that is the mean anomaly of the disturbed planet, and 
f m» are the mean anomalies of the disturbing planets m, j. 

In order to present more conveniently the equations defining the various 
parts of Sx^, Sy^, Sz^, let us put 


V aaf/ffOc + !l0 

Z-f 8a ) \ 9p 2x 9(3' 

a » P J X 


G 2y ) dt 


1' 


0 


9x o 9y o 

9 L 0 J » 9 io“ C 2 x ‘ — ° 2 y) dt ' J l 9 L^ ” 2 x ' 9 L 7 v, 2 y 


- |S-fl 


9x 


G„ + 


9y r 


dt 


%tj\w a a 2* + ^; %' dt 


ia_ |Y 3x o 

■*<oJW u 2y 


9y c 

G„ G_ I dt 


(21) 


for q = x, y. With this abbreviated notation, we have 


6c l2mjk (mjk) 


-z k/[:-? 

a , (3 


»Yr 


5T G 2xmjk (mjk) + ST G 2ymjk (mjk) 


dt 


n//[ 


9x 9y 

U r* / -\ ■ w 


3 n" -SM I I 8^ G 2xm jk (mjk) + 8^ G 2ymjk (m j k) 


dt 2 , 


where q = x, y, and 


,2 2mjk (mjk) * % J«1 ° 


'/ C 


2zmjk (m j k) dt - V «Z G 2zmjk (m j k) dt * (22) 
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The integrands on the right side of equations (22) can now be developed in 

Fourier series. Integrating these Fourier representations, we obtain other 

Fourier representations of the integrals. Multiplying these Fourier 

representations by the Fourier series representing the coefficients 3q/3a, 

2 -4 

-3 p Lq (Sq/Sco^), q^, -q^, and adding the different results, we obtain the 

Fourier representation 6x^ m j^(mjk), ^ Z 2mjk^ m ^^' We note that 

the constant coefficient in the Fourier series representing these different 

integrands, i. e. , the coefficients of the argument 0, will, when integrated 

once, give rise to a numerical coefficient multiplied by t; when integrated 

2 

twice, it will give rise to a numerical coefficient multiplied by t . Hence, 

the final representations of Sq^^-^Cmjk) will contain, besides the purely 

periodic terms given by the Fourier representations in three arguments, 

mixed terms composed of the time t multiplied by Fourier series in one 

2 

argument and, in the case of q = x, y only, the square of the time (t ) 
multiplied by Fourier series in one argument. These mixed terms will be 
added to the perturbations ^^kt^* ^2ktt^' # 


A computer program has been constructed with the series ^xmjk^'M* 
G 2 ymjk( m fk)> ^ 2 zmjk^ m ^k) as i n P u ^ and, as output, Fourier representations 
^ X 2mjk^ m ^^ Sy 2 m jk(mjk), ^ z 2mjk^ m ^^ an< 3 corre .spending mixed terms 

in 6x 2kt (k), &y 2kt (k), 6z 2kt (k), 6x 2ktt (k), Sy^k). 


For the evaluation of 5 
recall the following relations: 



jYf dt = t j f dt - j*J f dt 2 , 

J J t f dt 2 = t JJf dt 2 - 2 J/J f dt 3 , 


x, y, z, we must 


(23) 


where f is any function of time t. 
given by 


The equations defining Sq^^Cjk) will be 
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Sq 2jk ( J k ) = 


2>ffe 

ti da J L 9(3 


9y ( 

G z x jk^ k) + 3jT 


G 2yjk (jk) 


dt 


* 2 j 

3 p L. 


-4 8q ff[ 9x 0 
0 9 “oJJL 9w o 


By 




0 


Sxjk'M + 3^ G 2yj k« k > 


dt 



9y C 

G 2xjkt (:ik) + Bp" 


G 2yjkt (jk) 


dt 


+ 2 


where q = x, y, and 



■ 9y o 

G_ -1 *(jk) + G ? •! * 
2xjkt J 3 cOq 2yjkt 


(jk) 


dt“ 


6z 2jk (jk) = q 2 J q l °2zjk ( j k) dt “ q l J q 2 G 2zjk (jk) dt 

" q 2 JJ q l G 2 Z jkt (jk)dt2+ q l JJ <*2 G 2zjkt (jk) dfc2 • 


(24) 


Through the multiplication- of- series approach or the double-harmonic- 

analysis technique, we can develop the Fourier representations of all the 

integrands appearing in the above equations and then evaluate the Fourier 

representations of 8q^ j^jk) for q= x, y, z. We note again that the constant 

terms in the various harmonic representations of the above integrands will 

2 

give rise to mixed terms with coefficients t and t in the expressions for 

3 

8q 2 j k (jk), for q = x, y, z. Mixed terms with coefficient t will also appear in 
the cases for q = x, y because of the presence of triple integrals. These 
various mixed terms appearing in Sq^j^jk) will be included in the perturba^- 
tions 6 q 2kt (k), oq 2ktt (k)> and 6< l2kttt^ for q= x » y ’ z ' 
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The equations defining Sq^j^tjk) will be given by 


q 2jkt ( j k) = Z) 

a, P 


9x o . 9 y 0 

0(3 G 2xjkt^ k) + 9ff" G 2yjkt^ k) 


dt 


2 -4 ^9- 

3 (i L q 3» 0 


JP 


C 2xjkt (jk) + 


»y, 


G ,„ iM (:i k > 


9co Q 2yjkt 


dt 


where q = x, y, and 


6z 2jkt (jk) = q 2 | q l G 2zjkt (jk) dt " q l J q 2 G 2 Z jkt (jk) dt 


(25) 


Again, through the double -harmonic- analysis technique or the multiplication- 

of-series approach, we can get the harmonic representations of 6q« (jk) 

^jkt 

for q = x, y, z. Also, we expect mixed terms with coefficient t in Sq (jk) 

2 Zjkt 

for q = x, y, z and with coefficient t in the case of q = x, y. Since 6 q^ jR^jk) 

is already multiplied by t, these mixed terms will have coefficients 
2 3 

t and t . As before, these mixed terms will be included in the perturbations 
6< >2ktt (k) - 5< >2kttt (k) - 


Finally, for the evaluation of Sc^kt^’ 5q 2ktt^ k ^ ^ q 2kttt^ for q= x ’ V' z 
and ^^ktttt^ ^ or q = Xj Y> we raust recall the following relations: 

J*t 2 f dt = t 2 Jf dt - 2 tjjf dt 2 + 2jJJf dt 3 , 


JJ t 2 f dt = t 2 JJf dt 2 - 4t JJJf dt 3 + 6 jjjy f dt^ 


(26) 


where f is any function of time t. 
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We mentioned earlier the contributions to 6 ' 9 q 2kttt^ 

obtained while we were deriving expressions for 5 ^mjk^ 1 * 1 ^’ Sq^j^jk), 
6 q 2 jkt(j k )‘ a< ^ition t° tkese contributions, we have the following: 


6 W k) = - 2 


S“JJ[ 


9x o 9y o 

8(3“ G 2 xktt (k) + 8 fT G 2 yktt v 


i 


(k) dt‘ 


+ 4 3 (i L 


2 -4 8 q_ 


0 8oj 


;)I 


9x 0 9y 0 

G_ . /k) + -r — — G_ . „ 

8 oJq 2 xktt 8 (j 0 q 2 yktt 


dt~ 


where q = x, y; 


6z 2kt (k) = - 2 «2 IJll G 2zktt (k) dt2 + 2 q l JJ % G 2zktt (k) 


dt' , (27) 


6 q 2 ktt (k) = S la ( [ 

a, (3 J L 


9x 0 9 y 0 

8 p~ G 2 xktt (k) + 0 p“ G 2 yktt (k) 


dt 




3 11 ^0 G 2xktt <k| + 3u 0 G 2yktt* k) J 


dt 2 , 


where q = x, y; and 


6 z 2 ktt* k) ~ q 2 J q l G 2 zktt^ k) dt " q l | q 2 G 2 zktt* k) dt 


(28) 


3 4 

We note that mixed terms with coefficients t and t will appear from the 
expressions of ^^kt^’ dq 2 ktt^ §* ven ec L ua tions (27) and (28). These terms 
can be added to those defining S q 2 j c ^(k), ^^ktttt^' Terms that are purely 
periodic and are expressed in Fourier series in one argument will appear 
and are given by bq^^k), where 
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5q 2k M 


“ESffitS 


G 2xktt (k > + ST G 2 y ktt (k) 


dt' 


- 6 

for q = x, y, and 
6z 2k (k) = 2q 2 




3 “ 2 l ° 4 %)M b 


G 2xktt (k) + 8 co q G 2yktt (k ^ 


dt 4 , 


G 2^t« dt3 - 2 ^ 


G 2zktt< k > dt ■ 


(29) 


3 4 

Again, the terms may give rise to mixed terms with coefficients t and t . 
These will be added to ^^kttt^ ^ or q = X * ^ ’ Z anc * *° ^ q 2ktttt^ ^ or q = x> Y’ 


A computer program THEORY 2 has been constructed to compute 

5q 2jk^ k ^ 5q 2jkt^ k ^’ 5q 2k^’ 6q 2kt^* * * * * 5q 2ktttt^' The m P ut of this 
program is j, k. The final output is the Fourier representations of these 

perturbations. In this program, we followed the double- and single -harmonic - 

analysis methods; we did not apply the multiplication- of- series technique. 
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5. NUMERICAL APPLICATION 


In the previous section, we outlined the method followed for computing 
the second-order perturbations in 5x, 6y, and Sz. We have two main com- 
puter programs. By use of the harmonic -analysis approach, program 
THEORY 2 computes the periodic and secular perturbations expressed in 
Fourier series in the two mean anomalies L, f ^ of the disturbing and the 
disturbed planets and also in one mean anomaly of the disturbed planet. 

The second main program computes the periodic perturbations expressed 

in Fourier series in the three mean anomalies: l , S. ., the mean anomalies 

m’ y 

of the disturbing planets, and the mean anomaly of the disturbed planet. 

We used the multiplication-of-series approach, which required carrying 
the multiplication to a certain tolerance. This tolerance is taken to be di- 
rectly proportional to the divisor when we compute the integrand that will 
be integrated once. For the case of the integrand that will be integrated 

twice, we take the tolerance to be directly proportional to the square of the 

- 1 3 

divisor (the constant of proportionality is 10 ). This variable tolerance 

device will assure us that there has been no loss of any significant digits 
owing to the small divisor. 

The author will soon publish the details of these two main programs and 
the different subroutines associated with them. 


31 




6. NUMERICAL RESULTS 


In this section, we present the results of the computation of the second- 
order perturbation of Mars containing the masses of Jupiter and Saturn; 
according to the notation given previously, we give the results of ^^mjk 
(q = x, y, z), where m =6, j = 5, and k - 4. The other results will be given 
in another paper. 


Tables 1, 2, and 3 give the periodic part of the Fourier series repre- 
sentation of 5x~ .. , 8y , and dz . The mixed terms arising from 

Zmjk 7 2mjk 2mjk 

the evaluation of these perturbations are given in Tables 4, 5, and 6. The 

-1 9 

coefficients in these mixed terms are computed up to a tolerance of 10 
As a check, it would be interesting to compare the series we obtained for 
these perturbations with the results obtained from numerical integration. 


The above series representations (periodic and mixed) are simply the 
analytical solution of the set of differential equations (2), where G ^ G ^ , 


and G 


2 z 


are replaced by G 0 > G 0 , and G 0 . 

2xmjk 2ymjk Zzmjk 


That set of differ- 


ential equations has been solved numerically, using Cowell's method of 
numerical integration. In applying this method, the tenth difference has been 
neglected and the interval of integrations is taken to be 10 days. The initial 
values of the numerical integration are chosen such that 5x, 6y, 5 z at 
t = 90 days and t = 100 days are given by the analytical solutions of 5x, 5y, 6 z. 


The evaluation of G^xmjk* ^2ymjk’ ^2 zmjk numer i ca l integration 

of the differential equation was carried out by use of the original definition 


of these G's as given by equations (17). The integration has been carried 
out up to t = 40, 000 days. 
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When we compare the results of the numerical integration with the 
analytical representation, deviation is found between the two. The deviations 
found in the comparison of the perturbations in x and y are periodic in char- 

acter, with the amplitude increasing with time. The amplitude reaches 

_7 

5X10 around t = 20, 000. The deviation found in the comparison of the 
perturbation in z is again periodic, with smaller amplitude. The amplitude 

-9 

reaches 1 X 10 . The disagreement between the numerical solution and the 

analytical representation of the perturbation in x, y is very alarming. 

However, we must expect a satisfactory agreement if the starting values 

used in initiating the numerical integration are given to a great accuracy. 

These starting values have been obtained, as mentioned earlier, from the 

analytical solution. In obtaining the analytical solution, we carried the 

-13 

evaluation of the different integrals involved up to a tolerance of 10 ; i.e., 

- 1 3 

terms with absolute values less than 10 have been neglected. These 

terms may add up, causing the accuracy of the evaluation of the integral to 
-13 

be more than 10 . We must remember, also, that these integrals must be 

multiplied by the partial derivatives 9 xq/9oJq, ,8y /0.<o , 9 xq/9.Lq, 9Yq/9Lq, . . 

The coefficients of the harmonic representation of these partial derivatives 
2 

amount to 10 . Thus, the accuracy of the evaluation of the periodic repre- 
sentation of the perturbation in x, y, and z may amount to 10 ^ or even 


Numerical integration of the differential equations defining the second- 
order perturbation is very sensitive to the starting values, which we have 
just found may be in error to within 10 to 10 . To meet that situation, 

we can apply differential corrections to the starting values, such that the 
deviation between the numerical integration and the analytical solution is 
minimum, in the least- squares sense. 


The author owes this idea to Prof. G. M. Clemence. 
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We have applied differential corrections where our equation of condition 
corresponded to deviations at t = 2 00, 1800, 3400, . . . , 19,400 days. The 
results of this follow: 

At t = 90 days, 

A (6x) = -8. 6198240974 X 10" U 
A (6y) = -1. 3774726475 X 10" 10 
A (6z) = +3. 2748239038 X 10 -10 . 

At t = 1 00 days, 

A (6x) = -3. 6212757640 X 10 -11 
A (Sy) = -2. 4084833490 X lO' 10 
A (6z) = +3. 5234609123 X 10" 10 . 

When we apply these corrections to the starting values of the numerical 

integration, the agreement between the analytical solution and the numerical 

integration improves appreciably. The deviation, after the integration is 

-10 -9 

carried to 40, 000 days, never exceeds 4 X10 in x, y and 1x10 in z, 
an excellent agreement indeed. This comparison is shown in Table 7. 
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Table 1. Fourier representation of ^ x 2 m jk (P erio< ^ ic P art )- The coefficients are in units of 10 
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Table 2. Fourier representation of (periodic part). The coefficients are in units of 10 
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Table 3. Fourier representation of ^ z 2 m j^ (periodic part). The coefficients are in units of 10 
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g ir^oio n* *n ft » ho»m l un I • m t • -*in» oj t 

■g it m fV *4 • Nr. t | | ft •» 


«on«m«Ar(r 

O I t I lN»«w 
O I ^ 

t 


N«NOr*nm«o»(H»* 4 *<«o* 4 * 4 «> 4 ^>«N*«oiAnhoo- 
-<NNN Nin «0 4 >r< 14 ■« ft ft Mftmri-««|r 4 r 4 

* * •* rt i • • > ft ft 


^ OHNWN^OHlMNHO»tNlft^«N^OHffN^O- 4 NHO*<OrtinftffNH( 

lift I ft ft I ft ft i lift ft ft ft ft III- 


i t ft ft ft • ft ft ft ft 


aoN»i» 4 NnKiM 90 ffH#coo(c^ cm^Q 4 , rgmo «'fO<nn 
« n >o I Mno • im<nm ft meorsJl^r | biownh i r« 4 pj 43 -h | eg 4 rg 4- l i -« 

M 14 " l-H I N | I I r* ft I || | • N H ft 

<« I II 


«f ^ lO « ^ N dlCMfl - 
» ft) m i m »-* r- —• 
o I I I tn 


^N«wo>t 
-4- mm i 

— » r-t i 

ft 


N 4 NNOH 

I m rg »4 h 

I r-» 


rg 1 ft fM m (C O>rt «* 

I I 

I in | 
ft 


43 m o l I I 4 h-ON I 4 

N r< 4 g 4 I 

ft I 


_« 0 -gNWrt 0 HN<*ig^ 0 HNmgNr* 0 r*N»ngNH 0 HNHig<fiN^ 0 HNmgmNtg 0 HNmgm 

l l ft ft ii ill ill i» 

o'O'O'O'tooo®® ccoor-t'-r-r-r-r- > o««^' 0 >o<(}iri^inif\ia^^g gggggggmmmmmmmNN 
l l l l l l l l l l i l 11 l l l t till l » l l l t l I I I 1 I III 1 l I I l l l l l t 

^ 43 43 vfl 43 g g) g) iO 43 'O 43 43 41 4 ) 43 4 ) 4 ) 43 43 43 43 43 <0 43 43 43 43 41 43 43 43 43 43 43 43 4 ! 4 4 ) 43 4 4 ) 4 ) 43 4 43 4 


a rj h h h o r> oo h in 
c m rg | } m h n ini (? 4 
■“ I i i inoom g i 

« i I i i 


NO 43 CX 3 r~ ■— • O 
I N in go n I r* 
-4 43 >1- I 

-♦ I in 


O'r»h.n®a 4 ommo 4 9 'oonmginoo h 4 Nh fmn 
l ini r* o m rg m 4 ec r- in i rg m o- m 4 r-* 1 iwiigm 
r-i -g 4 1 iNineo 1 43 m 1 1 Nn 1 

|hn 1 g h 1 1 

1 


rg m 
ffl <0 m 
o 
o 


o oownnn «Ninmo «mmr4r-f^4>g43«f»>43gMmr-r*«>m>rHrgcMOP“43g«3(Mocom o 
M*imOH r- r- 43 43 1 a>Ni<H»iNNO 1 KiniMCC i t 4- m ac cq .g rg m — 

r-l 43 rg O' 4N I NHlfiOUlH I miftHtn | j OH |pg 

H I rg H 04 >gn 1 I 4 H | I 

I IW N I | 


mgi 3 m(MH 0 niMmgm{M- 40 'gNmgi<iNH 0 r*Ni*iggmNH 0 -iNingmNn 0 r 4 ing«iNH 0 

* ill III lilt lilt « I I » I » ft 1 1 1 

m mmg 4 - g 4 4 - gggrtKimmKimmNNNNNNNNiMHHHHH-g^oooooooHrtHHHw 

m l l l I 1 1 1 1 1 1 1 1 l 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

v£> 4 ^ -4 4 ‘ 4 ’ gggggggg ggggggggggggggggggggggggggggggggggg 


-4 g con 
® m <ng 
in in n 
rg IT. CM 

I 


O 

■4 


4 »emM\i 
cnoinn n g 
I lOO n O 
I Oh'l 


4 ) 


f* mONrlH 

O' O' r- m m 
m in cr m 
Hcfl 4 4 


in OISONNONO-O O' eg 
m in I ri4MnriNNN 
II Hh aw n«l I 

rg m ao | 

4 - I 


n n O' MnniNin^ H 

I n | NO Oinn 

-* r- rg n 
I in 

l 


g 

i 1 in 


>r o f‘ 0 'inNin 0 ' 0 ' 0 <inNmNO 0 'N 0 '®inrMn 4 hgNHinH 0 ' 0 ' 4 ingo«gogNN®^cingH 

a g Nonoo>N o> 43 4 a- m O' m no «3 in g 4 r- cm 4 om 43 mo— omroo— *— < »— • t — • o 
o 4 n o 4 00 oinoo'in 4 O' —• cc in cm O' in 43 O — 1 H-<gNin 1 I -* 

O t <— 1 O' t •— « cm 4 43 4 rtmm I gnn 1 1 rg -g 1 I 

rg pr, rg rg in I IT I | 

« » « 


HO^Ncng nnN-*Og Nn 4 - in 4 m rg -g O «-g rg K'gingmNgonNm 4 gmNgo-Nm 4 in 4 m 

• III till till till till 


ir> 4 - ggggggm(nm(n«uf|i«immNNNNNNNNN(MgggH* 4 gnggoOOOOOOO l -in -g rg 

• ft I 1 ft I I ft I I ft I I I « I I I I I ft 1 ft I I ft ft I ft I ft j ft I I ft 


gjNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN NN <M rg 


71 



Table 4. Fourier representation of Sx^ (mixed part) 


T measured from epoch in days. 

-19 

The coefficients are in units of 10 
The argument is 1 ^ mean anomaly of Mars* 
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sin 
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-26 T 2 
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-78999480 T 
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-8840577 T 

3 

-14383 T 

-977485 T 

4 

-1786 T 

-110676 T 

5 

-217 T 

-12768 T 

6 

-26 T 

-1493 T 

7 

-3 T 

-176 T 

8 

0 T 

-21 T 

9 

0 T 

-3 T 
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Table 5, Fourier representation of Sy^ (mixed part) 


T measured from epoch in days. 

The coefficients are in units of 1 0 * 

The argument is kf mean anomaly of Mars. 
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0 T 
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3 T 
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Table 6. Fourier representation of Sz^ (mixed part) 

T measured from epoch in days. 

_ 1 9 

The coefficients are in units of 10 
The argument is kf i ^ mean anomaly of Mars. 
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-1 T 

0 T 
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Table 7. Comparison of numerical integration with the analytical solution. 
The values are expressed in units of 10^ 
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